% Calibration
clc
clear
close all

addpath functions

% SN: Nordhaus view, constant tax rate
% CN: Nordhaus view, increasing tax rate
% SW: Disaster view, no mean reversion
% CW: Disaster view, mean reversion
% CWQ: Disaster view, mean reversion, qualitative

caslist={'CW','SW','CWQ','SN','CN'};
for jjc=1:length(caslist)
    
    cas=caslist{jjc};
    
    clearvars -except cas jjc caslist
    
    if strcmp(cas,'SN')
        
        mu=0.02;
        rho=0;
        delta=0.99;
        phi=0;
        
        xi=0.21;
        gamma=10;
        
        eta=3;
        psi=0;
        omega=0;
        
        lambdabar=0.03;
        alpha=0;
        chi=0;
        v=0;
        
        mulambda=lambdabar*(1-alpha-chi*xi);
        
        % Horizon length of term structure
        maxmat=1000;
        mm1=1000;
        mm2=50;
        
        mud=mu+(eta-1)*lambdabar*xi;
        mux=-lambdabar*phi*xi;
        muy=-lambdabar*psi*xi;
        
        
        % Damages
        muz=mud;
        piz=0;
        etaz=-eta/10;
        
    elseif strcmp(cas,'SW')
        
        mu=0.02;
        rho=0;
        delta=0.99;
        phi=0;
        
        xi=0.21;
        gamma=10;
        
        eta=3;
        psi=0;
        omega=0;
        
        lambdabar=0.03;
        alpha=0;
        chi=0;
        v=0.1;
        
        mulambda=lambdabar*(1-alpha-chi*xi);
        
        % Horizon length of term structure
        maxmat=1000;
        mm1=1000;
        mm2=50;
        
        mud=mu+(eta-1)*lambdabar*xi;
        mux=-lambdabar*phi*xi;
        muy=-lambdabar*psi*xi;
        
        
        % Damages
        muz=0;
        piz=0;
        etaz=eta/10;
        
    elseif strcmp(cas,'CN')
        
        mu=0.02;
        rho=0.85;
        delta=0.99;
        phi=0.025;
        
        xi=0.21;
        gamma=10;
        
        eta=3;
        psi=0.24;
        omega=.915;
        
        lambdabar=0.03;
        alpha=.75;
        chi=0.05;
        v=0.1;
        
        mulambda=lambdabar*(1-alpha-chi*xi);
        
        % Horizon length of term structure
        maxmat=1000;
        mm1=1000;
        mm2=50;
        
        mud=mu+(eta-1)*lambdabar*xi;
        mux=-lambdabar*phi*xi;
        muy=-lambdabar*psi*xi;
        
        
        % Damages
        muz=mud;
        piz=-1/10;
        etaz=-eta/10;
        
    elseif strcmp(cas,'CW')
        
        mu=0.02;
        rho=0.85;
        delta=0.99;
        phi=0.025;
        
        xi=0.21;
        gamma=10;
        
        eta=3;
        psi=0.24;
        omega=.915;
        
        lambdabar=0.03;
        alpha=.75;
        chi=0.05;
        v=0.1;
        
        mulambda=lambdabar*(1-alpha-chi*xi);
        
        % Horizon length of term structure
        maxmat=1000;
        mm1=1000;
        mm2=50;
        
        mud=mu+(eta-1)*lambdabar*xi;
        mux=-lambdabar*phi*xi;
        muy=-lambdabar*psi*xi;
        
        
        % Damages
        muz=[mud-2*lambdabar*xi*eta mud-2*lambdabar*xi*eta mud-2*lambdabar*xi*eta];
        piz=[-1 -1/2 -1/10];
        etaz=[eta eta/2 eta/10];
        
    elseif strcmp(cas,'CWQ')
        
        mu=0.02;
        rho=0.85;
        delta=0.99;
        phi=0.025;
        
        xi=0.21;
        gamma=10;
        
        eta=3;
        psi=0.24;
        omega=.915;
        
        lambdabar=0.03;
        alpha=.75;
        chi=0.05;
        v=0.1;
        
        mulambda=lambdabar*(1-alpha-chi*xi);
        
        % Horizon length of term structure
        maxmat=1000;
        mm1=1000;
        mm2=50;
       
        mud=mu+(eta-1)*lambdabar*xi;
        mux=-lambdabar*phi*xi;
        muy=-lambdabar*psi*xi;
        
        % Damages
        muz=0;
        piz=-1/10;
        etaz=eta/10; % (for dividends: eta)
       
    end
    
    
    
    % How many calibrations
    ndamages=length(muz);
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Paths of the economy: examples
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % No disaster
    dc_t(1,1)=mu-lambdabar*xi;
    dd_t(1,1)=mud-eta*lambdabar*xi;
    lambda_t(1:10,1)=lambdabar;
    x_t(1,1)=0;
    y_t(1,1)=0;
    
    % Convergence to steady state with no disasters
    for t=2:1000
        dc_t(t,1)=mu+x_t(t-1);
        x_t(t,1)=mux+rho*x_t(t-1);
        dd_t(t,1)=mud+y_t(t-1);
        y_t(t,1)=muy+omega*y_t(t-1);
        lambda_t(t,1)=mulambda+alpha*lambda_t(t-1)+v*x_t(t-1);
    end
    
    dclr=dc_t(1000,1);
    xlr=x_t(1000,1);
    ylr=y_t(1000,1);
    ddlr=dd_t(1000,1);
    lambdalr=lambda_t(1000,1);
    
    % No disaster case, but increase in lambda_t through change in xt and
    % yt
    dc_t(1000,1)=dclr;
    x_t(1000,1)=xlr;
    y_t(1000,1)=ylr;
    dd_t(1000,1)=ddlr;
    lambda_t(1000,1)=lambdalr;
    x_t(1000)=x_t(1000)+0.03;
    y_t(1000)=y_t(1000)+0.03;
    for t=1001:2000
        dc_t(t,1)=mu+x_t(t-1);
        x_t(t,1)=mux+rho*x_t(t-1);
        dd_t(t,1)=mud+y_t(t-1);
        y_t(t,1)=muy+omega*y_t(t-1);
        lambda_t(t,1)=mulambda+alpha*lambda_t(t-1)+v*x_t(t-1);
    end
    
    if strcmp(cas,'CW')
        figure
        subplot(2,1,1);
        plot(zeros(61,1),'k--','LineWidth',2);
        hold on
        plot(log(lambda_t(990:1050)/lambda_t(999)),'k','LineWidth',2);
        legend('Trend Growth','Above-Trend Growth','location','NorthEast')
        title('Disaster probability \lambda_{t}: log deviation from long-run mean')
        xlim([1 60])
        xlabel('Years')
        ylabel('\lambda_{t} Log Deviations')
        subplot(2,1,2);
        plot(ddlr*(1:61),'k--','LineWidth',2)
        hold on
        plot(cumsum(dd_t(990:1050)),'k','LineWidth',2);
        legend('Trend Growth','Above-Trend Growth','location','SouthEast')
        title('Rents')
        xlabel('Years')
        xlim([1 60])
        ylabel('Rents')
        set(gcf, 'units', 'inches', 'position', [2 2 12 7])
        set(gcf, 'Color', 'w');
        save2pdf('figures/path_nodisaster');
        close
    end
    
    ll_old=lambda_t;
    
    dd_t_nojump=dd_t;
    
    % Increase in xt and yt -- then disaster -- then recovery
    dc_t(1000,1)=dclr;
    dd_t(1000,1)=ddlr;
    x_t(1000,1)=xlr;
    y_t(1000,1)=ylr;
    lambda_t(1000,1)=lambdalr;
    x_t(1000)=x_t(1000)+0.03;
    y_t(1000)=y_t(1000)+0.03;
    for t=1001:2000
        if t~=1015
            dc_t(t,1)=mu+x_t(t-1);
            x_t(t,1)=mux+rho*x_t(t-1);
            dd_t(t,1)=mud+y_t(t-1);
            y_t(t,1)=muy+omega*y_t(t-1);
            lambda_t(t,1)=mulambda+alpha*lambda_t(t-1)+v*x_t(t-1);
        elseif t==1015 % Disaster
            dc_t(t,1)=mu+x_t(t-1)-xi;
            x_t(t,1)=mux+rho*x_t(t-1)+xi*phi;
            dd_t(t,1)=mud+y_t(t-1)-eta*xi;
            y_t(t,1)=muy+omega*y_t(t-1)+xi*psi;
            lambda_t(t,1)=mulambda+alpha*lambda_t(t-1)+v*x_t(t-1)+chi*xi;
        end
    end
    
    
    if strcmp(cas,'CW')
        figure
        subplot(2,1,1);
        plot(log(ll_old(990:1050)/ll_old(999)),'k--','LineWidth',2);
        hold on
        plot(log(lambda_t(990:1050)/lambda_t(999)),'k','LineWidth',2);
        legend('Above-Trend Growth / No Disaster','Above-Trend Growth / With Disaster','location','NorthEast')
        title('Disaster Probability \lambda_{t}: Log Deviation From Long-Run Mean')
        xlim([1 60])
        xlabel('Years')
        ylabel('\lambda_{t} Log Deviations')
        subplot(2,1,2);
        hold on
        plot(cumsum(dd_t_nojump(990:1050)),'k--','LineWidth',2);
        plot(cumsum(dd_t(990:1050)),'k','LineWidth',2);
        xlim([1 60])
        legend('Above-Trend Growth / No Disaster','Above-Trend Growth / With Disaster','location','SouthEast')
        xlabel('Years')
        ylabel('Rents')
        title('Rents')
        set(gcf, 'units', 'inches', 'position', [2 2 12 7])
        set(gcf, 'Color', 'w');
        save2pdf('figures/path_disaster');
        close
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Model solution
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Consumption strips
    a=nan(maxmat,1);
    b=nan(maxmat,1);
    f=nan(maxmat,1);
    a(1)=log(delta) + (1-gamma)*mu + log( (1-lambdabar) + lambdabar*exp(-(1-gamma)*xi));
    b(1)=1-gamma;
    f(1)=(exp(-(1-gamma)*xi)-1)/( 1 + lambdabar*(exp(-(1-gamma)*xi)-1));
    
    for j=2:maxmat
        a(j)=log(delta)+(1-gamma)*mu + a(j-1) + b(j-1)*mux + f(j-1)*(mulambda + (alpha-1)*lambdabar) ...
            + log((1-lambdabar) + lambdabar*exp((-(1-gamma)+b(j-1)*phi + f(j-1)*chi)*xi));
        b(j)=(1-gamma)+b(j-1)*rho + f(j-1)*v;
        f(j)=f(j-1)*alpha + (exp((-(1-gamma)+b(j-1)*phi + f(j-1)*chi)*xi) - 1) / (1 + lambdabar*(exp((-(1-gamma)+b(j-1)*phi ...
            + f(j-1)*chi)*xi)-1));
    end
    
    % Dividend strips
    ad=nan(maxmat,1);
    bd=nan(maxmat,1);
    ed=nan(maxmat,1);
    fd=nan(maxmat,1);
    
    ad(1)=log(delta)-gamma*mu + mud + log((1-lambdabar) + lambdabar*exp((gamma-eta)*xi));
    bd(1)=-gamma;
    ed(1)=1;
    fd(1)=(exp((gamma-eta)*xi)-1)/(1+lambdabar*(exp((gamma-eta)*xi)-1));
    
    for j=2:maxmat
        ad(j)=log(delta)-gamma*mu + mud + ad(j-1) + bd(j-1)*mux + ed(j-1)*muy + fd(j-1)*(mulambda+lambdabar*(alpha-1)) ...
            +log((1-lambdabar) + lambdabar*exp((gamma-eta+bd(j-1)*phi+ed(j-1)*psi+fd(j-1)*chi)*xi));
        bd(j)=-gamma+bd(j-1)*rho+fd(j-1)*v;
        ed(j)=1+ed(j-1)*omega;
        fd(j)=alpha*fd(j-1)+ (exp((gamma-eta+bd(j-1)*phi+ed(j-1)*psi+fd(j-1)*chi)*xi)-1)/(1+lambdabar*(exp((gamma-eta+bd(j-1)*phi+ed(j-1)*psi+fd(j-1)*chi)*xi)-1));
    end
    
    % Damage strips
    az=nan(maxmat,ndamages);
    bz=nan(maxmat,ndamages);
    ez=nan(maxmat,ndamages);
    fz=nan(maxmat,ndamages);
    
    for jj=1:ndamages
        
        az(1,jj)=log(delta)-gamma*mu + muz(jj) + log((1-lambdabar) + lambdabar*exp((gamma+etaz(jj))*xi));
        bz(1,jj)=-gamma;
        ez(1,jj)=piz(jj);
        fz(1,jj)=(exp((gamma+etaz(jj))*xi)-1)/(1+lambdabar*(exp((gamma+etaz(jj))*xi)-1));
        
        for j=2:maxmat
            az(j,jj)=log(delta)-gamma*mu + muz(jj) + az(j-1,jj) + bz(j-1,jj)*mux + ez(j-1,jj)*muy + fz(j-1,jj)*(mulambda+lambdabar*(alpha-1)) ...
                +log((1-lambdabar) + lambdabar*exp((gamma+etaz(jj)+bz(j-1,jj)*phi+ez(j-1,jj)*psi+fz(j-1,jj)*chi)*xi));
            bz(j,jj)=-gamma+bz(j-1,jj)*rho+fz(j-1,jj)*v;
            ez(j,jj)=piz(jj)+ez(j-1,jj)*omega;
            fz(j,jj)=alpha*fz(j-1,jj)+ (exp((gamma+etaz(jj)+bz(j-1,jj)*phi+ez(j-1,jj)*psi+fz(j-1,jj)*chi)*xi)-1)/(1+lambdabar*(exp((gamma+etaz(jj)+bz(j-1,jj)*phi+ez(j-1,jj)*psi+fz(j-1,jj)*chi)*xi)-1));
        end
        
    end
    
    % Risk-free rate
    a_f=nan(maxmat,1);
    b_f=nan(maxmat,1);
    f_f=nan(maxmat,1);
    a_f(1)=log(delta) -gamma*mu + log( 1-lambdabar + lambdabar*exp(gamma*xi));
    b_f(1)=-gamma;
    f_f(1)=(exp(gamma*xi)-1)/(1+(exp(gamma*xi)-1)*lambdabar);
    
    for j=2:maxmat
        a_f(j)=log(delta)-gamma*mu + a_f(j-1) + b_f(j-1)*mux + f_f(j-1)*mulambda + (alpha-1)*f_f(j-1)*lambdabar ...
            + log(1-lambdabar + lambdabar*exp((gamma+b_f(j-1)*phi+f_f(j-1)*chi)*xi));
        b_f(j)=-gamma+b_f(j-1)*rho+f_f(j-1)*v;
        f_f(j)=f_f(j-1)*alpha + (exp((gamma+b_f(j-1)*phi+f_f(j-1)*chi)*xi)-1)/(1+lambdabar*(exp((gamma+b_f(j-1)*phi+f_f(j-1)*chi)*xi)-1));
    end
    
    
    %Implied discount rates, general case
    
    ngridpoints=200;
    
    % Consumption
    
    An_series=An(lambdabar,0,mm1,xi,chi,mulambda,alpha,v,mux,rho,phi,ngridpoints);
    
    logret=(1:mm1)'*mu + (mux*cumsum((1-rho.^(0:mm1-1))/(1-rho)))' + log(An_series) - a(1:mm1);
    
    % Dividends
    
    And_series=And(lambdabar,0,mm1,xi,chi,mulambda,alpha,v,mux,rho,phi,psi,omega,eta,ngridpoints);
    
    logretd=(1:mm1)'*mud + (muy*cumsum((1-omega.^(0:mm1-1))/(1-omega)))'+log(And_series) - ad(1:mm1);
    
    % Damages
    
    for jj=1:ndamages
        Anz_series(:,jj)=Anz(lambdabar,0,mm1,xi,chi,mulambda,alpha,v,mux,rho,phi,psi,omega,etaz(jj),piz(jj),ngridpoints);
        
        logretz(:,jj)=(1:mm1)'*muz(jj) + (piz(jj)*muy*cumsum((1-omega.^(0:mm1-1))/(1-omega)))'+log(Anz_series(:,jj)) - az(1:mm1,jj);
    end
    
    % Risk-free rate
    
    logretf = -a_f(1:mm1);
    

    % Implied leasehold discounts
    % ----------------------------------------
    
    pd_strips=exp(ad);
    pd_avg=cumsum(pd_strips);
    
    % Log price-dividend ratio of freehold
    logpf=log(pd_avg(end));
    
    % Discounts, UK bins
    clear logpd;
    logpd(1)=mean(log(pd_avg(80:99)));
    logpd(2)=mean(log(pd_avg(100:124)));
    logpd(3)=mean(log(pd_avg(125:149)));
    logpd(4)=mean(log(pd_avg(150:300)));
    logpd(5)=mean(log(pd_avg(701:999)));
    
    % Leasehold discounts, UK
    discounts_uk=logpd-logpf;
    
    % Discounts, SG bins
    clear logpd;
    logpd(1)=mean(log(pd_avg(50:70)));
    logpd(2)=mean(log(pd_avg(71:85)));
    logpd(3)=mean(log(pd_avg(86:90)));
    logpd(4)=mean(log(pd_avg(91:95)));
    logpd(5)=mean(log(pd_avg(96:100)));
    logpd(6)=mean(log(pd_avg(801:999)));
    
    % Leasehold discounts, SG
    discounts_sg=logpd-logpf;
    
    data_uk=[-0.176 -0.110 -0.089 -0.033 -0.003];
    data_sg=[-0.409 -0.275 -0.215 -0.148 -0.125 -0.01];
    labels_uk={'80-99','100-124','125-149','150-300','700+'};
    labels_sg={'50-70','71-85','86-90','91-95','96-100','800+'};
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Semi-Elasticity of PD Ratios of Different Properties to Disaster Risk
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    mult=0.01;
    eta_H=eta*(1+mult);
    eta_L=eta;
    mud_H=mu+(eta_H-1)*lambdabar*xi;
    mud_L=mu+(eta_L-1)*lambdabar*xi;
    
    % Dividend strips
    ad_H=nan(maxmat,1);
    bd_H=nan(maxmat,1);
    ed_H=nan(maxmat,1);
    fd_H=nan(maxmat,1);
    
    ad_L=nan(maxmat,1);
    bd_L=nan(maxmat,1);
    ed_L=nan(maxmat,1);
    fd_L=nan(maxmat,1);
    
    ad_H(1)=log(delta)-gamma*mu + mud_H + log((1-lambdabar) + lambdabar*exp((gamma-eta_H)*xi));
    bd_H(1)=-gamma;
    ed_H(1)=1;
    fd_H(1)=(exp((gamma-eta_H)*xi)-1)/(1+lambdabar*(exp((gamma-eta_H)*xi)-1));
    
    ad_L(1)=log(delta)-gamma*mu + mud_L + log((1-lambdabar) + lambdabar*exp((gamma-eta_L)*xi));
    bd_L(1)=-gamma;
    ed_L(1)=1;
    fd_L(1)=(exp((gamma-eta_L)*xi)-1)/(1+lambdabar*(exp((gamma-eta_L)*xi)-1));
    
    for j=2:maxmat
        ad_H(j)=log(delta)-gamma*mu + mud_H + ad_H(j-1) + bd_H(j-1)*mux + ed_H(j-1)*muy + fd_H(j-1)*(mulambda+lambdabar*(alpha-1)) ...
            +log((1-lambdabar) + lambdabar*exp((gamma-eta_H+bd_H(j-1)*phi+ed_H(j-1)*psi+fd_H(j-1)*chi)*xi));
        bd_H(j)=-gamma+bd_H(j-1)*rho+fd_H(j-1)*v;
        ed_H(j)=1+ed_H(j-1)*omega;
        fd_H(j)=alpha*fd_H(j-1)+ (exp((gamma-eta_H+bd_H(j-1)*phi+ed_H(j-1)*psi+fd_H(j-1)*chi)*xi)-1)/(1+lambdabar*(exp((gamma-eta_H+bd_H(j-1)*phi+ed_H(j-1)*psi+fd_H(j-1)*chi)*xi)-1));
        
        ad_L(j)=log(delta)-gamma*mu + mud_L + ad_L(j-1) + bd_L(j-1)*mux + ed_L(j-1)*muy + fd_L(j-1)*(mulambda+lambdabar*(alpha-1)) ...
            +log((1-lambdabar) + lambdabar*exp((gamma-eta_L+bd_L(j-1)*phi+ed_L(j-1)*psi+fd_L(j-1)*chi)*xi));
        bd_L(j)=-gamma+bd_L(j-1)*rho+fd_L(j-1)*v;
        ed_L(j)=1+ed_L(j-1)*omega;
        fd_L(j)=alpha*fd_L(j-1)+ (exp((gamma-eta_L+bd_L(j-1)*phi+ed_L(j-1)*psi+fd_L(j-1)*chi)*xi)-1)/(1+lambdabar*(exp((gamma-eta_L+bd_L(j-1)*phi+ed_L(j-1)*psi+fd_L(j-1)*chi)*xi)-1));
    end
    
    pd_strips_H=exp(ad_H);
    pd_avg_H=cumsum(pd_strips_H);
    freeholdpd_H=pd_avg_H(end);
    weights_H=pd_strips_H/freeholdpd_H;
    
    pd_strips_L=exp(ad_L);
    pd_avg_L=cumsum(pd_strips_L);
    freeholdpd_L=pd_avg_L(end);
    weights_L=pd_strips_L/freeholdpd_L;
    
    elast_H=weights_H'*fd_H;
    elast_L=weights_L'*fd_L;
    elast=elast_H-elast_L;
    
    if strcmp(cas,'CW')
        pd_avg_H(end)
        pd_avg_L(end)
        elast
    end
    
    % Conditional volatility of the bond return at different maturities
    volrb(1)=0; % Note that bf and ff for n=0 are both zero
    for j=2:maxmat
        volrb(j,1)=(b_f(j-1)*phi+f_f(j-1)*chi)^2 * lambdabar * (1-lambdabar) * xi^2;
    end
    % This is the volatility of the 1-year return
    sqrt(volrb(1:100));
    
    if strcmp(cas,'CW')
        
        figure
        hold on
        plot(logretf./(1:mm1)','r--','LineWidth',3)
        plot(logretd./(1:mm1)','b','LineWidth',2)
        for jj=1:ndamages
            if jj==1
                plot(logretz(1:mm1,3)./(1:mm1)','k','LineWidth',2)
            elseif jj==2
                plot(logretz(1:mm1,2)./(1:mm1)','k--','LineWidth',2)
            elseif jj==3
                plot(logretz(1:mm1,1)./(1:mm1)','k:','LineWidth',2)
            end
        end
        xlim([1 700])
        xlabel('Years')
        legend('Risk-Free','Housing','Damages (\theta=1/10)','Damages (\theta=1/2)','Damages (\theta=1)')
        axis tight
        ylabel('Discount Rate')
        set(gcf, 'units', 'inches', 'position', [2 2 12 7])
        set(gcf, 'Color', 'w');
        save2pdf('figures/rfdr');
        close
        
        figure
        colormap('jet')
        set(gca,'fontname','Times New Roman','fontSize',14);
        bar([discounts_uk;data_uk]')
        xlabel('Leasehold Maturity');
        ylabel('Discount');
        legend('Model','Data','Location','SouthEast')
        set(gca,'Xticklabel',labels_uk)

        title('UK');
        set(gcf, 'units', 'inches', 'position', [2 2 12 7])
        set(gcf, 'Color', 'w');
        save2pdf('figures/discounts_uk');
        close
        
        figure
        colormap('jet')
        set(gca,'fontname','Times New Roman','fontSize',14);
        bar([discounts_sg;data_sg]')
        xlabel('Leasehold Maturity');
        ylabel('Discount');
        legend('Model','Data','Location','SouthEast')
        set(gca,'Xticklabel',labels_sg)
        title('SG')
        set(gcf, 'units', 'inches', 'position', [2 2 12 7])
        set(gcf, 'Color', 'w');
        save2pdf('figures/discounts_sg');
        close
    end
    
    rpmat=(logretz-logretf)./(1:mm1)';
    save(['temp/' cas '.mat'],'rpmat');
    
    
end

% Comparison of different models

clear
figure
hold on

load 'temp/SN.mat'
plot(ones(size(rpmat))*0.015,'b--','LineWidth',2);
plot(rpmat+0.015,'k','LineWidth',2)

load 'temp/CN.mat'
plot(rpmat+0.015,'k--','LineWidth',2)

load 'temp/SW.mat'
plot(rpmat+0.015,'r','LineWidth',2)

load 'temp/CWQ.mat'
plot(rpmat+0.015,'r--','LineWidth',2)

xlim([1 700])
xlabel('Years')
ylim([-0.005 0.045])

ylabel('Discount Rate')

legend('Risk-Free Rate','Tax View of Climate Change (Constant)','Tax View of Climate Change (Increasing)','Disaster View of Climate Change (No Adaptation)','Disaster View of Climate Change (Adaptation)','location','southeast')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('figures/qual');
close

